*do file reproduces all figures from the online appendix

*************************************************************
*Figure A.1: Horizontal Well Production: Millions of $ of BOE
*************************************************************
use "$tmp/final_playonly.dta",clear
* RPI Index (SD Shift) regression
quietly reghdfe allprod_hwell ddsdrpiscore $rhs $main
lincom ddsdrpiscore
loc bstd = r(estimate)
loc sestd = r(se)
loc pstd = r(p)
loc stars ""
if `pstd' <= .1 {
loc stars "*"
}
if `pstd' <= .05 {
loc stars "**"
}
if `pstd' <= .01 {
loc stars "***"
}
loc stdlab = "RPI Index (SD Shift)× Post = `bstd'`stars' [`sestd']"
reghdfe allprod_hwell eyibin1_12_tvsdrpiscore_ts1 eyibin1_13_tvsdrpiscore_ts1 eyibin1_14_tvsdrpiscore_ts1 eyibin1_15_tvsdrpiscore_ts1 eyibin1_16_tvsdrpiscore_ts1 eyibin1_17_tvsdrpiscore_ts1 eyibin1_18_tvsdrpiscore_ts1 eyibin1_19_tvsdrpiscore_ts1 eyibin1_20_tvsdrpiscore_ts1 eyibin1_21_tvsdrpiscore_ts1 eyibin1_22_tvsdrpiscore_ts1 eyibin1_23_tvsdrpiscore_ts1 eyibin1_24_tvsdrpiscore_ts1 eyibin1_25_tvsdrpiscore_ts1 eyibin1_26_tvsdrpiscore_ts1 eyibin1_27_tvsdrpiscore_ts1 eyibin1_28_tvsdrpiscore_ts1 eyibin1_29_tvsdrpiscore_ts1 $rhs $main
estimates store rpi_results

* Above median Indicator regression  

quietly reghdfe allprod_hwell ddamedian1 $rhs $main
lincom ddamedian1
loc bMed = r(estimate)
loc seMed = r(se)
loc pMed = r(p)
loc stars ""
if `pMed' <= .1 {
loc stars "*"
}
if `pMed' <= .05 {
loc stars "**"
}
if `pMed' <= .01 {
loc stars "***"
}
loc Medlab = "Above Median IndicatorXPost = `bMed'`stars' [`seMed']"

reghdfe allprod_hwell eyibin1_12_tvamedian1_ts1 eyibin1_13_tvamedian1_ts1 eyibin1_14_tvamedian1_ts1 eyibin1_15_tvamedian1_ts1 eyibin1_16_tvamedian1_ts1 eyibin1_17_tvamedian1_ts1 eyibin1_18_tvamedian1_ts1 eyibin1_19_tvamedian1_ts1 eyibin1_20_tvamedian1_ts1 eyibin1_21_tvamedian1_ts1 eyibin1_22_tvamedian1_ts1 eyibin1_23_tvamedian1_ts1 eyibin1_24_tvamedian1_ts1 eyibin1_25_tvamedian1_ts1 eyibin1_26_tvamedian1_ts1 eyibin1_27_tvamedian1_ts1 eyibin1_28_tvamedian1_ts1 eyibin1_29_tvamedian1_ts1 $rhs $main
estimates store quartile_results


coefplot (rpi_results, keep(eyibin1_*_tvsdrpiscore_ts1) ///
          rename(eyibin1_12_tvsdrpiscore_ts1 = "-11" ///
                 eyibin1_13_tvsdrpiscore_ts1 = "-10" ///
                 eyibin1_14_tvsdrpiscore_ts1 = "-9" ///
                 eyibin1_15_tvsdrpiscore_ts1 = "-8" ///
                 eyibin1_16_tvsdrpiscore_ts1 = "-7" ///
                 eyibin1_17_tvsdrpiscore_ts1 = "-6" ///
                 eyibin1_18_tvsdrpiscore_ts1 = "-5" ///
                 eyibin1_19_tvsdrpiscore_ts1 = "-4" ///
                 eyibin1_20_tvsdrpiscore_ts1 = "-3" ///
                 eyibin1_21_tvsdrpiscore_ts1 = "-2" ///
                 eyibin1_22_tvsdrpiscore_ts1 = "-1" ///
                 eyibin1_23_tvsdrpiscore_ts1 = "0" ///
                 eyibin1_24_tvsdrpiscore_ts1 = "1" ///
                 eyibin1_25_tvsdrpiscore_ts1 = "2" ///
                 eyibin1_26_tvsdrpiscore_ts1 = "3" ///
                 eyibin1_27_tvsdrpiscore_ts1 = "4" ///
                 eyibin1_28_tvsdrpiscore_ts1 = "5" ///
                 eyibin1_29_tvsdrpiscore_ts1 = "6") ///
          color(orange) msymbol(circle) msize(medium)) ///
         (quartile_results, keep(eyibin1_*_tvamedian1_ts1) ///
          rename(eyibin1_12_tvamedian1_ts1 = "-11" ///
                 eyibin1_13_tvamedian1_ts1 = "-10" ///
                 eyibin1_14_tvamedian1_ts1 = "-9" ///
                 eyibin1_15_tvamedian1_ts1 = "-8" ///
                 eyibin1_16_tvamedian1_ts1 = "-7" ///
                 eyibin1_17_tvamedian1_ts1 = "-6" ///
                 eyibin1_18_tvamedian1_ts1 = "-5" ///
                 eyibin1_19_tvamedian1_ts1 = "-4" ///
                 eyibin1_20_tvamedian1_ts1 = "-3" ///
                 eyibin1_21_tvamedian1_ts1 = "-2" ///
                 eyibin1_22_tvamedian1_ts1 = "-1" ///
                 eyibin1_23_tvamedian1_ts1 = "0" ///
                 eyibin1_24_tvamedian1_ts1 = "1" ///
                 eyibin1_25_tvamedian1_ts1 = "2" ///
                 eyibin1_26_tvamedian1_ts1 = "3" ///
                 eyibin1_27_tvamedian1_ts1 = "4" ///
                 eyibin1_28_tvamedian1_ts1 = "5" ///
                 eyibin1_29_tvamedian1_ts1 = "6") ///
          color(navy) msymbol(diamond) msize(medium)), ///
    vertical yline(0, lcolor(black) lpattern(dash)) ///
    xline(0, lcolor(black) lpattern(solid)) ///
    ytitle("Coefficient Value") ///
    xtitle("Year Relative to Fracking's Introduction") ///
    title("Horizontal Well Production: Millions of $ of BOE") ///
    legend(label(2 "`stdlab'") ///
           label(4 "`Medlab'") ///
           pos(6) cols(1)) ///
    scheme(s1mono) graphregion(color(white))

graph export "$outappend/figA1.pdf", replace 


*************************************************************
*Figure A2: Age-Adjusted Overall Mortality per 100,000
**Variable naming in code (see readme for more details)
**adr2000dc0_gen0_eth00, Age-Adjusted overall death rate (dc0), 2000 population, both genders (gen0)
**adr2000dc0_gen1_eth00, Age-Adjusted overall death rate (dc0), 2000 population, male (gen1)
**adr2000dc0_gen2_eth00, Age-Adjusted overall death rate (dc0), 2000 population, female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval i = 0/2 {
	qui reghdfe adr2000dc0_gen`i'_eth00 ddtquartile1 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)

	reghdfe adr2000dc0_gen`i'_eth00 alpha_* $rhs  $main 
	eventstudy gen`i' 0 `beta' `se' `p' 
	// Add panel letters: A=overall, B=male, C=female  
	local panel = cond(`i'==0, "A", cond(`i'==1, "B", "C"))
	graph export "$outappend/figA2`panel'.pdf", replace 
}


*************************************************************
*Figure A3: Crude Overall Mortality Rate per 100,000: No Controls
**Requires Restricted CDC Mortality Data
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen1_eth00_a2564, crude overall death rate (dc0), male (gen1)
**cdrdc0_gen2_eth00_a2564, crude overall death rate (dc0), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval i = 0/2 {
	qui reghdfe cdrdc0_gen`i'_eth00_a2564 ddtquartile1 $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)

	reghdfe cdrdc0_gen`i'_eth00_a2564 alpha_* $main 
	eventstudy gen`i' 0 `beta' `se' `p' 
	// Add panel letters: A=overall, B=male, C=female  
	local panel = cond(`i'==0, "A", cond(`i'==1, "B", "C"))
	graph export "$outappend/figA3`panel'.pdf", replace 
}


*************************************************************
*Figure A5: Overall Mortality Effects by Age
**Requires Restricted CDC Mortality Data
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a0024, crude overall death rate (dc0), both genders (gen0), under 25 (a0024)
**cdrdc0_gen0_eth00_a2544, crude overall death rate (dc0), both genders (gen0), ages 25-44 (a2544)
**cdrdc0_gen0_eth00_a4564, crude overall death rate (dc0),both genders (gen0), ages 45-64 (a4564)
**cdrdc0_gen0_eth00_a6599, crude overall death rate (dc0),both genders (gen0), ages 65+ (a6599)

*Figure A6: Overall Mortality Effects by Age: Male
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a0024, crude overall death rate (dc0), male (gen1), under 25 (a0024)
**cdrdc0_gen0_eth00_a2544, crude overall death rate (dc0), male (gen1), ages 25-44 (a2544)
**cdrdc0_gen0_eth00_a4564, crude overall death rate (dc0), male (gen1), ages 45-64 (a4564)
**cdrdc0_gen0_eth00_a6599, crude overall death rate (dc0), male (gen1), ages 65+ (a6599)

*Figure A7: Overall Mortality Effects by Age: Female
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a0024, crude overall death rate (dc0), female (gen2), under 25 (a0024)
**cdrdc0_gen0_eth00_a2544, crude overall death rate (dc0), female (gen2), ages 25-44 (a2544)
**cdrdc0_gen0_eth00_a4564, crude overall death rate (dc0), female (gen2), ages 45-64 (a4564)
**cdrdc0_gen0_eth00_a6599, crude overall death rate (dc0), female (gen2), ages 65+ (a6599)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval i = 0/2{
	foreach a in 0024 2544 4564 6599 {
		reghdfe cdrdc0_gen`i'_eth00_a`a' ddtquartile1 $rhs $main
		lincom ddtquartile1
		loc beta: di %9.3fc r(estimate)
		loc se: di %9.3fc r(se)
		loc p = r(p)
		reghdfe cdrdc0_gen`i'_eth00_a`a'  alpha_* $rhs  $main 
		eventstudy age`i'`a' 0 `beta' `se' `p'
		loc figure = cond(`i'==0,"A5",cond(`i'==1,"A6","A7"))
		loc panel = cond("`a'"=="0024","A",cond("`a'"=="2544","B",cond("`a'"=="4564","C","D")))
		graph export "$outappend/fig`figure'`panel'.pdf", as(pdf) name("Graph") replace
}
}

*************************************************************
*Figure A8: Crude Death Rates: Internal and External Causes of Death by Gender
**Variable naming in code (See readme for more details)
**cdrdc30_gen1_eth00_a2564, crude internal death rate (dc30),male (gen1)
**cdrdc30_gen1_eth00_a2564, crude internal death rate (dc30), male (gen1)
**cdrd31_gen2_eth00_a2564, crude external death rate (d31), female (gen2)
**cdrdc31_gen2_eth00_a2564, crude external death rate (d31), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

foreach d in 30 31 {
	forval i = 1/2{
		reghdfe cdrdc`d'_gen`i'_eth00_a2564 ddtquartile1 $rhs $main
		lincom ddtquartile1
		loc beta: di %9.3fc r(estimate)
		loc se: di %9.3fc r(se)
		loc p = r(p)
		reghdfe cdrdc`d'_gen`i'_eth00_a2564  alpha_* $rhs  $main 
		*set panels and ranges
		if `d' == 30 {
			loc panel = cond(`i' == 1,"A","B")
			loc ymin = -40
			loc yint = 20
			loc ymax = 40		
		}
		else if `d' == 31 {
			loc panel = cond(`i' == 1,"C","D")
			loc ymin = -20
			loc yint = 10
			loc ymax = 20	
		}
		eventstudy dc`d'gen`i' 0 `beta' `se' `p' o `ymin' `yint' `ymax'
		graph export "$outappend/figA8`panel'.pdf", as(pdf) name("Graph") replace
	}
}

*************************************************************
*Figure A9: External Causes of Death: Deaths of Despair
**Variable naming in code: 
**cdrdc1_gen0_eth00_a2564, crude suicide death rate (dc1), both genders (gen0)
**cdrdc2_gen0_eth00_a2564, crude alcohol related death rate (dc2), both genders (gen0)
**cdrdc3_gen0_eth00_a2564, crude drug overdose related death rate (dc3), both genders (gen0)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval d = 1/3 {
	qui reghdfe cdrdc`d'_gen0_eth00_a2564 ddtquartile1 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)
	reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* $rhs  $main 
	eventstudy dc`d' 0 `beta' `se' `p' o -12 6 12
	local panel = cond(`d' == 3,"A",cond(`d' == 1,"B","C"))
	graph export "$outappend/figA9`panel'.pdf", replace 
}

*************************************************************
*Poisson Regressions: 
*Figure A.4: Overall Mortality
**Variable naming in code (See readme for more details)
**Panel A: ndc0_gen0_eth00_a2564, Raw death counts, overall mortality (dc0), both genders (gen0)
**Panel B: ndc0_gen1_eth00_a2564, Raw death counts, overall mortality  (dc0), male (gen1)
**Panel C: ndc0_gen2_eth00_a2564, Raw death counts, overall mortality  (dc0), female (gen2)

*Figure A.10: Internal vs. External Causes of Death
**Panel A: ndc30_gen0_eth00_a2564, Raw death counts, internal mortality (dc30), both genders (gen0)
**Panel B: ndc31_gen0_eth00_a2564, Raw death counts, external mortality (dc31), both genders (gen0)

*Figure A.11: Internal and External Causes of Death by Gender

**Panel A: ndc30_gen1_eth00_a2564,Raw death counts, internal mortality (dc30), male (gen1)
**Panel B: ndc30_gen2_eth00_a2564, Raw death counts, internal mortality (dc30), female (gen2)
**Panel C: ndc31_gen1_eth00_a2564,Raw death counts, external mortality (dc31), male (gen1)
**Panel D: ndc31_gen2_eth00_a2564, Raw death counts, external mortality (dc31), female (gen2)

*Figure A.12: Cardiovascular vs Other Internal Causes of Death
**Panel A: ndc4_gen0_eth00_a2564, Raw death counts, internal mortality (dc30), both genders (gen0)
**Panel B: ndc32_gen0_eth00_a2564, Raw death counts, external mortality (dc31), both genders (gen0)

*Figure A.13: Deaths of Despair
**Panel A: ndc3_gen0_eth00_a2564, Raw death counts, drug overdose mortality (dc3), both genders (gen0)
**Panel B: ndc1_gen1_eth00_a2564, Raw death counts, suicide mortality (dc1), both genders (gen0)
**Panel C: ndc2_gen2_eth00_a2564, Raw death counts, alcohol related mortality (dc2), both genders (gen0)

*************************************************************
use "$tmp/final_playonly.dta",clear

foreach d in 0 1 2 3 4 30 31 32 {
	forval i = 0/2 {
		//Define figure numbers and panel letters based on death type (d) and gender (i)
		if `d' == 0 {  //Overall mortality - Figure A.4
			local figure "A4"
			if `i' == 0 local panel "A"
			if `i' == 1 local panel "B" 
			if `i' == 2 local panel "C"
		}
		else if `d' == 30 & `i' == 0 {  //Internal mortality, both genders - Figure A.10, Panel A
			local figure "A10"
			local panel "A"
		}
		else if `d' == 31 & `i' == 0 {  //External mortality, both genders - Figure A.10, Panel B
			local figure "A10"
			local panel "B"
		}
		else if `d' == 30 & `i' == 1 {  //Internal mortality, male - Figure A.11, Panel A
			local figure "A11"
			local panel "A"
		}
		else if `d' == 30 & `i' == 2 {  //Internal mortality, female - Figure A.11, Panel B
			local figure "A11"
			local panel "B"
		}
		else if `d' == 31 & `i' == 1 {  //External mortality, male - Figure A.11, Panel C
			local figure "A11"
			local panel "C"
		}
		else if `d' == 31 & `i' == 2 {  //External mortality, female - Figure A.11, Panel D
			local figure "A11"
			local panel "D"
		}
		else if `d' == 4 & `i' == 0 {  //Cardiovascular mortality - Figure A.12, Panel A
			local figure "A12"
			local panel "A"
		}
		else if `d' == 32 & `i' == 0 {  //Non-cardiovascular internal mortality - Figure A.12, Panel B
			local figure "A12"
			local panel "B"
		}
		else if `d' == 3 & `i' == 0 {  //Drug overdose mortality - Figure A.13, Panel A
			local figure "A13"
			local panel "A"
		}
		else if `d' == 1 & `i' == 0 {  //Suicide mortality - Figure A.13, Panel B
			local figure "A13"
			local panel "B"
		}
		else if `d' == 2 & `i' == 0 {  //Alcohol mortality - Figure A.13, Panel C
			local figure "A13"
			local panel "C"
		}
		else {
			continue  //Skip combinations that don't correspond to any figure - no regression runs
		}

	
		ppmlhdfe ndc`d'_gen`i'_eth00_a2564 ddtquartile1 $rhs   if balancesample == 1, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000) exposure( pop_gen`i'_eth00_a2564 )
		nlcom exp(_b[ddtquartile1])-1
		mat b = r(b)
		mat V = r(V)
		loc beta_1: di %9.3f b[1,1]
		loc beta = trim("`beta_1'")
		loc se_1: di %4.3f sqrt(V[1,1])
		loc se = trim("`se_1'")
		loc p =  2*normal(-abs(b[1,1]/sqrt(V[1,1])))
		loc stars ""
		if `p' <= .1 {
			loc stars "*"
		}
		if `p' <= .05 {
			loc stars "**"
		}
		if `p' <= .01 {
			loc stars "***"
		}
	
	ppmlhdfe ndc`d'_gen`i'_eth00_a2564  alpha_* $rhs  if balancesample == 1, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000) exposure( pop_gen`i'_eth00_a2564 )
        // Initialize min/max values
        local min_val = .
        local max_val = .
        
        // Get coefficient names
        matrix b = e(b)
        local colnames : colnames b
        
        // Loop through alpha coefficients
        foreach v of local colnames {
            // Only process alpha_12 through alpha_29
            if regexm("`v'", "^alpha_(1[2-9]|2[0-9])$") {
                // Get coefficient and SE
                local coef = _b[`v']
                local s = _se[`v']
                
                // Calculate transformed CI bounds (exp(b±1.96*se)-1)
                local ci_low = exp(`coef' - 1.96*`s') - 1
                local ci_high = exp(`coef' + 1.96*`s') - 1
                
                // Update min/max
                if `min_val' == . {
                    local min_val = `ci_low'
                    local max_val = `ci_high'
                }
                else {
                    local min_val = min(`min_val', `ci_low')
                    local max_val = max(`max_val', `ci_high')
                }
            }
        }
		        local abs_max = max(abs(`min_val'), abs(`max_val'))
	
	 *Calculate symmetric bounds
     local abs_max = max(abs(`min_val'), abs(`max_val'))
        

	*Round to nice values for percentage scales
	if `abs_max' < 0.1 {
		local base = 0.05   // For very small effects (<10%)
	}
	else if `abs_max' < 0.25 {
		local base = 0.1    //For small effects (10-25%)
	}
	else if `abs_max' < 0.5 {
		local base = 0.25   // For medium effects (25-50%)
	}
	else if `abs_max' < 1 {
		local base = 0.5    // For larger effects (50-100%)
	}
	else if `abs_max' < 2 {
		local base = 1      // For effects between 100-200%
	}
	else if `abs_max' < 5 {
		local base = 2.5    // For effects between 200-500%
	}
	else {
		local base = 5      // For very large effects (>500%)
	}

	local y_max = ceil(`abs_max'/`base')*`base'

	* Choose interval (smaller intervals for percentage scales)
	if `y_max' > 5 {
		local y_int = 2.5
	}
	else if `y_max' > 2 {
		local y_int = 1
	}
	else if `y_max' > 1 {
		local y_int = 0.5
	}
	else if `y_max' > 0.5 {
		local y_int = 0.25
	}
	else if `y_max' > 0.25 {
		local y_int = 0.1
	}
		else {
			local y_int = 0.05
		}        
        * Create y-axis options for coefplot
        local yaxis_opts "ylabel(-`y_max'(`y_int')`y_max') yscale(range(-`y_max' `y_max'))"
		coefplot, transform(*=exp(@)-1) vertical keep(alpha*) omit order(alpha12-alpha29) graphregion(color(white)) lcolor(white) lwidth(vthick) ciopts(lcolor(black)) mcolor(black) msymbol(o) yline(0, lc(gs6) lp(dash)) xline(11, lc(black)) ytitle("Coefficient Estimate") subtitle("Top-Quartile × Post = `beta'`stars' [`se']", size(medsmall) color(black))     xtitle("Year Relative to Fracking's Introduction") xlabel(1 "-11" 2 "-10" 3 "-9" 4 "-8" 5 "-7" 6 "-6" 7 "-5" 8 "-4" 9 "-3" 10 "-2" 11 "-1" 12 "0" 13 "1" 14 "2" 15 "3" 16 "4" 17 "5" 18 "6")  ysc(titleg(2)) legend(off) `yaxis_opts'
		graph export "$outappend/fig`figure'`panel'.pdf", as(pdf) name("Graph") replace
}
}


*************************************************************
*Figure A14: Internal and External Causes of Death: Differences by Gender
*************************************************************

use "$tmp/final_playonly.dta",clear
*Panel A: Internal Causes of Death - Men vs Women 
* Define internal causes
local death_causes "4 5 6 7 8 9 10"
local death_labels `""cardiovascular deaths" "cancer deaths" "respiratory deaths" "infection-related deaths" "brain disease deaths" "kidney/urethra deaths" "nutrition-related deaths""'
* Create storage for results
tempname results1
postfile `results1' outcome_num gender_num coef se str50 outcome_label str20 gender_label mean_val n_obs using temp_results1, replace
local outcome_count = 0

* Loop through each internal cause of death
foreach dc_num in `death_causes' {
    local ++outcome_count    
    * Get outcome label 
    local outcome_lab : word `outcome_count' of `death_labels' 
   * Men (gen1) regression
    preserve
    local men_var "cdrdc`dc_num'_gen1_eth00_a2564"
    * Calculate baseline mean for men
    qui sum `men_var' if post == 1 & tquartile1 == 0
    local men_mean = r(mean)
    * Run regression for men
    qui reghdfe `men_var' ddtquartile1 $rhs $main
    post `results1' (`outcome_count') (1) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`outcome_lab'") ("Men") (`men_mean') (e(N))
    restore
    
    * Women (gen2) regression  
    preserve
    local women_var "cdrdc`dc_num'_gen2_eth00_a2564"
    
    * Calculate baseline mean for women
    qui sum `women_var' if post == 1 & tquartile1 == 0
    local women_mean = r(mean)
    * Run regression for women
    qui reghdfe `women_var' ddtquartile1 $rhs $main
    post `results1' (`outcome_count') (2) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`outcome_lab'") ("Women") (`women_mean') (e(N))
    restore
}

* Create Figure 
postclose `results1'
use temp_results1, clear

* Create confidence intervals
gen ci_hi = coef + 1.96 * se
gen ci_lo = coef - 1.96 * se

* Create y-axis positioning
gen snum = outcome_num
replace snum = outcome_num + 0.08 if gender_num == 1  // Men slightly above
replace snum = outcome_num - 0.08 if gender_num == 2  // Women slightly below

* Create labels with means
gen y_label = outcome_label + " [Mean = " + string(mean_val, "%9.0f") + "]" if gender_num == 1
replace y_label = "[Mean = " + string(mean_val, "%9.0f") + "]" if gender_num == 2

* Set up y-axis labels
local ylab ""
levelsof snum, local(slist)
foreach s in `slist' {
    preserve
    keep if abs(snum - `s') < 0.0001
    if _N > 0 {
        local lab = y_label[1]
        if "`lab'" != "" & "`lab'" != "." {
            local ylab "`ylab' `s' "`lab'""
        }
    }
    restore
}

* Create Panel A plot
twoway (rspike ci_lo ci_hi snum if gender_num == 1, horizontal lc(blue)) ///
       (scatter snum coef if gender_num == 1, mc(blue) ms(circle)) ///
       (rspike ci_lo ci_hi snum if gender_num == 2, horizontal lc(red)) ///
       (scatter snum coef if gender_num == 2, mc(red) ms(circle)), ///
       graphregion(color(white)) ///
       title("Panel A: Internal Causes", size(medium)) ///
       subtitle("Outcome: deaths per 100K (ages 25-64)", size(small)) ///
       ytitle("") ysc(titleg(3)) ///
       ylabel(`ylab', angle(0) labsize(small) nogrid) ///
       xtitle("Post × Top-Quartile Coefficient", size(medium)) ///
       xlabel(, grid labsize(small)) xsc(titleg(3)) ///
       legend(order(2 "Men" 4 "Women") position(6) rows(1) size(medium)) ///
       plotregion(margin(l+5 r+5 t+3 b+3))

graph export "$outappend\figA14A.pdf", as(pdf) name("Graph") replace

* Clean up 
erase temp_results1.dta

* Load original data for Panel B


//Panel B: External Causes of Death - Men vs Women 
use "$tmp/final_playonly.dta",clear
* Define external causes 
local death_causes "1 2 3 11 12 13"
local death_labels `""suicide deaths" "alcohol-related deaths" "drug-related deaths" "motor vehicle deaths" "other accidents deaths" "homicide deaths""'

* Create storage for results
tempname results2
postfile `results2' outcome_num gender_num coef se str50 outcome_label str20 gender_label mean_val n_obs using temp_results2, replace

local outcome_count = 0

* Loop through each external cause of death
foreach dc_num in `death_causes' {
    local ++outcome_count
    
    * Get outcome label 
    local outcome_lab : word `outcome_count' of `death_labels'
    
    * Men (gen1) regression
    preserve
    local men_var "cdrdc`dc_num'_gen1_eth00_a2564"
    
    * Calculate baseline mean for men
    qui sum `men_var' if post == 1 & tquartile1 == 0
    local men_mean = r(mean)
    
    * Run regression for men
    qui reghdfe `men_var' ddtquartile1 $rhs $main 
    post `results2' (`outcome_count') (1) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`outcome_lab'") ("Men") (`men_mean') (e(N))
    restore
    
    * Women (gen2) regression  
    preserve
    local women_var "cdrdc`dc_num'_gen2_eth00_a2564"
    
    * Calculate baseline mean for women
    qui sum `women_var' if post == 1 & tquartile1 == 0
    local women_mean = r(mean)
    
    * Run regression for women
    qui reghdfe `women_var' ddtquartile1 $rhs $main
    post `results2' (`outcome_count') (2) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`outcome_lab'") ("Women") (`women_mean') (e(N))
    restore
}

* Create Figure 
postclose `results2'
use temp_results2, clear

* Create confidence intervals
gen ci_hi = coef + 1.96 * se
gen ci_lo = coef - 1.96 * se

* Create y-axis positioning
gen snum = outcome_num
replace snum = outcome_num + 0.08 if gender_num == 1  // Men slightly above
replace snum = outcome_num - 0.08 if gender_num == 2  // Women slightly below

* Create labels with means
gen y_label = outcome_label + " [Mean = " + string(mean_val, "%9.0f") + "]" if gender_num == 1
replace y_label = "[Mean = " + string(mean_val, "%9.0f") + "]" if gender_num == 2

* Set up y-axis labels
local ylab ""
levelsof snum, local(slist)
foreach s in `slist' {
    preserve
    keep if abs(snum - `s') < 0.0001
    if _N > 0 {
        local lab = y_label[1]
        if "`lab'" != "" & "`lab'" != "." {
            local ylab "`ylab' `s' "`lab'""
        }
    }
    restore
}

* Create Panel B plot
twoway (rspike ci_lo ci_hi snum if gender_num == 1, horizontal lc(blue)) ///
       (scatter snum coef if gender_num == 1, mc(blue) ms(circle)) ///
       (rspike ci_lo ci_hi snum if gender_num == 2, horizontal lc(red)) ///
       (scatter snum coef if gender_num == 2, mc(red) ms(circle)), ///
       graphregion(color(white)) ///
       title("Panel B: External Causes", size(medium)) ///
       subtitle("Outcome: deaths per 100K (ages 25-64)", size(small)) ///
       ytitle("") ysc(titleg(3)) ///
       ylabel(`ylab', angle(0) labsize(small) nogrid) ///
       xtitle("Post × Top-Quartile Coefficient", size(medium)) ///
       xlabel(, grid labsize(small)) xsc(titleg(3)) ///
       legend(order(2 "Men" 4 "Women") position(6) rows(1) size(medium)) ///
       plotregion(margin(l+5 r+5 t+3 b+3))

graph export "$outappend\figA14B.pdf", as(pdf) name("Graph") replace
* Clean up
erase temp_results2.dta

*************************************************************
*Figure A15: SNAP Benefits
*************************************************************
use "$tmp\SNAP.dta",clear
reghdfe pct_snap ddtquartile1 $rhs, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)
	loc stars ""
	if `p' <= .1 {
		loc stars "*"
	}
	if `p' <= .05 {
		loc stars "**"
	}
	if `p' <= .01 {
		loc stars "***"
	}
	
reghdfe pct_snap  alpha_* $rhs , absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
coefplot, vertical keep(alpha*) omit order(alpha4-alpha11) graphregion(color(white)) lcolor(white) lwidth(vthick) ciopts(lcolor(black)) mcolor(black) msymbol(o) yline(0, lc(gs6) lp(dash)) xline(11, lc(black)) ytitle("Coefficient Estimate") subtitle("Top-Quartile × Post = `beta'`stars' [`se']", size(medsmall) color(black))     xtitle("Year Relative to Fracking's Introduction") xlabel(1 "-4" 2 "-3" 3 "-2" 4 "-1" 5 "0" 6 "1" 7 "2" 8 "3")  ysc(titleg(2) r(-.02 .02))  ylabel(-.02(.01).02) legend(off)
graph export "$outappend/figA15.pdf", as(pdf) name("Graph") replace


*************************************************************
*Figure A16: Overall Mortality: Controlling for Compositional Changesr
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen1_eth00_a2564, crude overall death rate (dc0), male (gen1)
**cdrdc0_gen2_eth00_a2564, crude overall death rate (dc0), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval i = 0/2 {
	reghdfe cdrdc0_gen`i'_eth00_a2564  ddtquartile1 share_gen`i'_eth00_a0024 share_gen`i'_eth00_a2544 share_gen`i'_eth00_a4564 share_gen`i'_eth00_a6599 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)

	reghdfe cdrdc0_gen`i'_eth00_a2564 alpha_* share_gen`i'_eth00_a0024 share_gen`i'_eth00_a2544 share_gen`i'_eth00_a4564 share_gen`i'_eth00_a6599 $rhs $main
	eventstudy gen`i' 0 `beta' `se' `p' o -40 10 40
	// Add panel letters: A=overall, B=male, C=female  
	local panel = cond(`i'==0, "A", cond(`i'==1, "B", "C"))
	graph export "$outappend/figA16`panel'.pdf", as(pdf) name("Graph") replace
}

*************************************************************
*Figure A17: Internal and External Causes of Death: Controlling for Compositional Changes
**Variable naming in code (See readme for more details)
**cdrdc30_gen1_eth00_a2564, crude internal death rate (dc30),male (gen1)
**cdrdc30_gen1_eth00_a2564, crude internal death rate (dc30), male (gen1)
**cdrd31_gen2_eth00_a2564, crude external death rate (d31), female (gen2)
**cdrdc31_gen2_eth00_a2564, crude external death rate (d31), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

foreach d in 30 31 {
	forval i = 1/2{
		reghdfe cdrdc`d'_gen`i'_eth00_a2564  ddtquartile1 share_gen`i'_eth00_a0024 share_gen`i'_eth00_a2544 share_gen`i'_eth00_a4564 share_gen`i'_eth00_a6599 $rhs $main
		lincom ddtquartile1
		loc beta: di %9.3fc r(estimate)
		loc se: di %9.3fc r(se)
		loc p = r(p)
		reghdfe cdrdc`d'_gen`i'_eth00_a2564 alpha_* share_gen`i'_eth00_a0024 share_gen`i'_eth00_a2544 share_gen`i'_eth00_a4564 share_gen`i'_eth00_a6599 $rhs $main
		*set panels and ranges
		if `d' == 30 {
			loc panel = cond(`i' == 1,"A","B")
			loc ymin = -40
			loc yint = 10
			loc ymax = 40		
		}
		else if `d' == 31 {
			loc panel = cond(`i' == 1,"C","D")
			loc ymin = -20
			loc yint = 5
			loc ymax = 20	
		}
		eventstudy dc`d'gen`i' 0 `beta' `se' `p' o `ymin' `yint' `ymax'
		graph export "$outappend/figA17`panel'.pdf", as(pdf) name("Graph") replace
	}
}

*************************************************************
*Figure A18: Overall Mortality by Gender: No Bakken Shale Play
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen1_eth00_a2564, crude overall death rate (dc0), male (gen1)
**cdrdc0_gen2_eth00_a2564, crude overall death rate (dc0), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

forval i = 0/2 {
	reghdfe cdrdc0_gen`i'_eth00_a2564  ddtquartile1 $rhs $restrict
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)

	reghdfe cdrdc0_gen`i'_eth00_a2564 alpha_*  $rhs $restrict
	eventstudy gen`i' 0 `beta' `se' `p' 
	// Add panel letters: A=overall, B=male, C=female  
	local panel = cond(`i'==0, "A", cond(`i'==1, "B", "C"))
	graph export "$outappend/figA18`panel'.pdf", as(pdf) name("Graph") replace
}

*************************************************************
*Figure A19: Internal vs. External Causes of Death: No Bakken Shale Play
**Variable naming in code: 
**cdrdc30_gen0_eth00_a2564, crude internal death rate (dc30), both genders (gen0)
**cdrdc31_gen1_eth00_a2564, crude external death rate (dc31), both genders (gen0)
*************************************************************
use "$tmp/final_playonly.dta",clear

*internal causes of death
foreach d in 30 31 {
	qui reghdfe cdrdc`d'_gen0_eth00_a2564 ddtquartile1 $rhs $restrict
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)
	reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* $rhs  $restrict 
	if `d' == 30 {
		loc panel "A"
		loc ymin = -40
		loc yint = 10
		loc ymax = 40		
	}
	else if `d' == 31 {
		loc panel "B"
		loc ymin = -20
		loc yint = 5
		loc ymax = 20				
	}
	eventstudy dc`d' 0 `beta' `se' `p' o `ymin' `yint' `ymax'
	graph export "$outappend/figA19`panel'.pdf", replace 
}


*************************************************************
*Figure A.20: Men/Women Working-Age Mortality Robustness
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
*************************************************************
use "$tmp/final_playonly.dta",clear
* Create storage for results using postfile
tempname results
postfile `results' ymn ysn festimate fse str50 flabym str50 flabys using temp_results, replace

//Sample 1: Balanced Sample, No ND or MT (FIRST)
preserve
keep if sfips2000 != "38" & sfips2000 != "30" & balancesample == 1
* Get sample mean for labeling
qui sum cdrdc0_gen0_eth00_a2564
local sample1_mean = r(mean)
local sample1_n = r(N)
local mean1_label = "Mean = " + string(`sample1_mean', "%9.2f")
local n1_label = "N = " + string(`sample1_n', "%9.0fc")
* Specification 1: No Controls, No Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (10) (1) (_b[ddtquartile1]) (_se[ddtquartile1]) ("Balanced Sample, No ND or MT") ("No Controls, No Weights")

* Specification 2: Controls, No Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (10) (2) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`mean1_label'") ("Controls, No Weights")

* Specification 3: No Controls, Weights  
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (10) (3) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`n1_label'") ("No Controls, Weights")

* Specification 4: Controls, Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (10) (4) (_b[ddtquartile1]) (_se[ddtquartile1]) ("") ("Controls, Weights")
restore

//Sample 2: Balanced Sample, No Restrictions (SECOND)
preserve
keep if balancesample == 1
* Get sample mean for labeling
qui sum cdrdc0_gen0_eth00_a2564
local sample2_mean = r(mean)
local sample2_n = r(N)
local mean2_label = "Mean = " + string(`sample2_mean', "%9.2f")
local n2_label = "N = " + string(`sample2_n', "%9.0fc")

* Specification 1: No Controls, No Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (20) (1) (_b[ddtquartile1]) (_se[ddtquartile1]) ("Balanced Sample, No Restrictions") ("No Controls, No Weights")

* Specification 2: Controls, No Weights  
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (20) (2) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`mean2_label'") ("Controls, No Weights")

* Specification 3: No Controls, Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (20) (3) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`n2_label'") ("No Controls, Weights")

* Specification 4: Controls, Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (20) (4) (_b[ddtquartile1]) (_se[ddtquartile1]) ("") ("Controls, Weights")
restore

//Sample 3: Unbalanced Sample, No ND or MT (THIRD)
preserve
keep if sfips2000 != "38" & sfips2000 != "30"
* Get sample mean for labeling
qui sum cdrdc0_gen0_eth00_a2564
local sample3_mean = r(mean)
local sample3_n = r(N)
local mean3_label = "Mean = " + string(`sample3_mean', "%9.2f")
local n3_label = "N = " + string(`sample3_n', "%9.0fc")

* Specification 1: No Controls, No Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (30) (1) (_b[ddtquartile1]) (_se[ddtquartile1]) ("Unbalanced Sample, No ND or MT") ("No Controls, No Weights")

* Specification 2: Controls, No Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs, absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (30) (2) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`mean3_label'") ("Controls, No Weights")

* Specification 3: No Controls, Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (30) (3) (_b[ddtquartile1]) (_se[ddtquartile1]) ("`n3_label'") ("No Controls, Weights")

* Specification 4: Controls, Weights
qui reghdfe cdrdc0_gen0_eth00_a2564 ddtquartile1 $rhs [aw=pop2000], absorb(splay1_year_fe cntyfips2000) vce(cl cntyfips2000)
post `results' (30) (4) (_b[ddtquartile1]) (_se[ddtquartile1]) ("") ("Controls, Weights")
restore

//Create the figure 
postclose `results'
use temp_results, clear
* Destring coefficient and standard error variables
destring festimate fse, replace force
* Create confidence intervals
gen ci95hi = festimate + 1.96 * fse
gen ci95lo = festimate - 1.96 * fse 
* Create y-axis position variable 
gen snum = .
replace snum = 34 - (ymn/10 - 1)*10 - (ysn - 1) if ymn == 10
replace snum = 24 - (ymn/10 - 2)*10 - (ysn - 1) if ymn == 20  
replace snum = 14 - (ymn/10 - 3)*10 - (ysn - 1) if ymn == 30
* Set up y-axis labels 
local ylab ""
levelsof ymn, local(ymlist)
foreach ym in `ymlist' {
    forval ys = 1/4 {
        preserve
        keep if ymn == `ym' & ysn == `ys'
        if _N > 0 {
            qui sum snum
            local s = r(mean)
            local lab = flabym[1]
            if "`lab'" != "" & "`lab'" != " " & "`lab'" != "." {
                local ylab "`ylab' `s' "`lab'""
            }
        }
        restore
    }
}

* Define colors and markers 
local colorlist "dknavy cranberry dkgreen dkorange"
local mlist "d o t s"
* Build graph components
local fgraph ""
local legend_order ""
* Add dummy series for legend
local i = 0
forval ys = 1/4 {
    local ++i
    local color: word `i' of `colorlist'
    local marker: word `i' of `mlist'
    local fgraph "`fgraph' (connected fse fse if _n == 10000, m(`marker') lc(`color') mc(`color'))"
    * Build legend
    preserve
    keep if ysn == `ys'
    if _N > 0 {
        local spec_lab = flabys[1]
        local legend_order "`legend_order' `i' "`spec_lab'""
    }
    restore
}

* Add actual data points
levelsof ymn, local(ymlist)
foreach ym in `ymlist' {
    local i = 0 
    forval ys = 1/4 {
        local ++i         
        local color: word `i' of `colorlist'
        local marker: word `i' of `mlist'
        local condition "ymn == `ym' & ysn == `ys'"
        
        local fgraph "`fgraph' (rspike ci95lo ci95hi snum if `condition', horizontal lc(`color')) (scatter snum festimate if `condition', mc(`color') m(`marker'))"
    }
}

* Create the final graph 
twoway `fgraph', graphregion(color(white)) ///
title("All-Cause Death Rate", size(medlarge) color(black)) ///
subtitle("Men and Women, ages 25-64", size(medium) color(black)) ///
ytitle("") ysc(titleg(3)) ///
ylab(`ylab', valuelabel angle(0) labsize(small) nogrid) ///
xtitle("Top-Quartile × Post", size(medium)) xlabel(, grid labsize(small)) xsc(titleg(3)) ///
legend(order(`legend_order') position(6) rows(1) region(lwidth(none)) size(small) symxsize(8)) ///
plotregion(margin(l+5 r+5 t+3 b+3))

graph export "$outappend\figA20.pdf", as(pdf) name("Graph") replace
* Clean up
erase temp_results.dta

*************************************************************
*Figure A21: Mortality: Geographic/Place Based Factors
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen1_eth00_a2564, crude overall death rate (dc0), male (gen1)
**cdrdc0_gen2_eth00_a2564, crude overall death rate (dc0), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta",clear

*create "pre" mortality variables
foreach d in 0 30 31 {
	cap gen adr`d'1990 =  adr2000dc`d'_gen0_eth00 if year == 1990
	cap egen adr2000dc`d'_y1990 = mean(adr`d'1990),by(cntyfips2000 )

	reghdfe cdrdc`d'_gen0_eth00_a2564  ddtquartile1 c.adr2000dc`d'_y1990#i.year l_pop_gen0_eth00_a2564 $rhs $main
	lincom ddtquartile1
	loc beta: di %9.3fc r(estimate)
	loc se: di %9.3fc r(se)
	loc p = r(p)

	reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* c.adr2000dc`d'_y1990#i.year l_pop_gen0_eth00_a2564 $rhs $main
	eventstudy dc`d' 0 `beta' `se' `p'
	local panel = cond(`d'==0, "A", cond(`d'==30, "B", "C"))
	graph export "$outappend/figA21`panel'.pdf", as(pdf) name("Graph") replace
}

*************************************************************
*Figure A22: Potential Confounding: Comparing Shale to No-Shale Counties within State
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen30_eth00_a2564, crude overall death rate (dc30)
**cdrdc0_gen31_eth00_a2564, crude overall death rate (d31)
*************************************************************
use "$tmp\allstates.dta",clear
foreach d in 0 30 31 {	
reghdfe cdrdc`d'_gen0_eth00_a2564 ddoverplay $rhs [aw = pop2000] if balance == 1 , absorb(state_year cntyfips2000) vce(cl cntyfips2000)
lincom ddoverplay
loc beta: di %9.3fc r(estimate)
loc se: di %9.3fc r(se)
loc p = r(p)
reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* $rhs [aw = pop2000] if balance == 1 , absorb(state_year cntyfips2000) vce(cl cntyfips2000)
 stateevent dc`d' 0 `beta' `se' `p'
 
loc panel = cond(`d'==0,"A",cond(`d'==30,"B","C"))
 graph export "$outappend/figA22`panel'.pdf", as(pdf) name("Graph") replace

 }
